A phase-field model by an Ising machine and its application to the phase-separation structure of a diblock polymer

A novel model to be applied to next-generation accelerators, Ising machines, is formulated on the basis of the phase-field model of the phase-separation structure of a diblock polymer. Recently, Ising machines including quantum annealing machines, attract overwhelming attention as a technology that opens up future possibilities. On the other hand, the phase-field model has demonstrated its high performance in material development, though it takes a long time to achieve equilibrium. Although the convergence time problem might be solved by the next-generation accelerators, no solution has been proposed. In this study, we show the calculation of the phase-separation structure of a diblock polymer as the equilibrium state using phase-field model by an actual Ising machine. The proposed new model brings remarkable acceleration in obtaining the phase-separation structure. Our model can be solved on a large-scale quantum annealing machine. The significant acceleration of the phase-field simulation by the quantum technique pushes the material development to the next stage.

Diblock polymers are a widely used material for structures 1 . The microstructure of a diblock polymer changes depending on the synthetic conditions. The process of the change in microstructure is spinodal decomposition, and various phase-separated structures, such as lamellar, cubic, hexagonal, and gyroid structures, can be obtained 2 . Since this material structure has a strong influence on mechanical properties 3 , it is very important to predict which structure will be obtained experimentally.
To date, several methods have been developed to analyze the equilibrium state of diblock polymers, such as the molecular dynamics 4 Monte Carlo method 5 , the self-consistent field theory 6 and the phase-field model 7,8 . The phase-field model is a continuum model that expresses the interface with a smooth function using a variable called the order parameter. This model has the advantage of being free from solving complicated boundary value problems. The governing equation can be derived in a relatively simple form. Hence, it is widely used to reproduce various material structures 1,[9][10][11][12][13][14][15][16] . The governing equation of the phase-field model is formulated by variation of the energy functional. However, conventional analysis sometimes takes a very long time to obtain an equilibrium state, which is a problem when a large-scale simulation is required. Hence, various schemes have been proposed to solve this problem [17][18][19] , which is the key to accelerating material development.
On the other hand, next-generation accelerators including quantum computers are steadily developing. In particular, quantum annealing machines [20][21][22][23] can search for the minimum value of the objective function at high speed. They are often applied to combinational optimization problems involving the objective functions of the target. If the objective function can be replaced with energy expressing the target phenomenon, quantum acceleration of finding of equilibrium state acquisition can be expected by a physics simulation [24][25][26][27] .
In this study, we propose an objective function for exploring the phase-separated structure of diblock polymers by annealing, with an eye toward the use of quantum annealing. On the basis of the obtained results, we confirm the tendency of the microstructure and the simulation performance by comparison with the phase-field simulation. The approach is based on the use of a global optimization metaheuristic algorithm, called simulated annealing, to directly minimize the Helmholtz free energy instead of minimizing it analytically and then solving www.nature.com/scientificreports/ the resulting nonlinear, partial differential equation, i.e., the governing equation of the phase-field model. In this paper, we illustrate the use of simulated annealing in the solution of the phase-field model by applying it to the formation of a microstructure in a diblock polymer. First, we review the results of the conventional analysis method of the phase-field model and propose the Ising model for solving the phase-field model in the following section. Subsequently, we present the calculation results of the newly defined quantity for the Ising model. Ising machines are used as a dedicated accelerator for solving Ising model, and from the simulation by the Ising machine, we find remarkable results, which indicate the effectiveness of the proposed approach. Finally, we conclude with remarks on the possibilities of a much larger-scale simulation and on the application of the annealing method to other problems related to continuum mechanics.
To verify the problem of spinodal decomposition for diblock polymers, we carry out conventional phase-field simulations. A phase diagram of the phase-separated structure in the diblock polymer is shown in Fig. 1 2 . The phase-field model is composed of the gradient energy, Flory-Huggins interaction energy 28,29 , and Ohta-Kawasaki long-distance interaction energy 30 . The fineness of the phase-separation structure is controlled by the magnitude of the Ohta-Kawasaki energy. The governing equation is derived through the Cahn-Hilliard equation 30 and discretized by the finite difference method following the conventional manner. The mobility M in the analysis is employed as M = 1.0 s −1 . The analysis area of 32 µm × 32 µm is discretized into 64 × 64 grids under the periodic boundary condition.
The initial distribution is given by random noise around f = 0.5 , which is shown in Fig. 2a(i). First, we conduct a numerical analysis with the constant of long-distance interaction terms c 0 equal to zero. Figure 2 a(ii) shows the equilibrium state obtained by phase-field analysis. Polymer A is shown in blue, and polymer B is shown in red, being separated into two large phases. These large phases are considered to minimize the surface energy. Next, we change the parameter c 0 in the long-distance term derived from the Ohta-Kawasaki energy 30 . We confirm that the incremental error at the final computational step becomes less than 1.2 × 10 −6 in all cases. The incremental error is defined as the product of the time increment per step and the maximum increment of the order parameter per step among all the grids. In Fig. 2b, the separation pattern, which is different from that in Fig. 2a, becomes finer with increasing parameter c 0 , which is considered to be induced by the effect of the longdistance energy term of the Ohta-Kawasaki energy. Additionally, the phase-separation structure in Fig. 2a(ii) seems cubic or hexagonal, as in Fig. 1, while that in Fig. 2b(iii) is lamellar-like, as in Fig. 1. The effect becomes more significant as the effect of the Ohta-Kawasaki energy increases. Both simulations take approximately 5960 s.
We performed a new phase-field simulation using an Ising machine which solves one of the formalized combinatorial optimization problems, i.e., quadratic unconstrained binary optimization (QUBO). QUBO uses a set of binary states. The objective function, called the Hamiltonian, is composed only of the quadratic polynomials of the states. One of the advantages of Ising machines is that once the problem is well formalized as a QUBOstyle Hamiltonian, Ising machines can be used to solve the problem without paying attention to the details of numerical calculations, such as the discretization error of finite time steps and the solver implementation method. www.nature.com/scientificreports/ Our proposed method formalizes the problem of spinodal decomposition for diblock polymers into the QUBO-style Hamiltonian. This Hamiltonian is composed of four terms that represent different physical behaviors of diblock polymers: the sum amount conservation term, the Flory-Huggins interaction energy term, the gradient energy term and the Ohta-Kawasaki long-distance interaction energy term. The strength of the impact on each of the four terms is controlled by the parameters α F , α I , α A and α OK , respectively, and the parameter f determines the fraction of monomer. The strength of the Ohta-Kawasaki energy term controls the fineness of the phase-separation structure. Note that Ising machines efficiently find the global minimum solution of a given Hamiltonian without becoming trapped at a local minimum solution. Thus, if a desired final state of a phase-field model is a local minimum free energy near the initial state, Ising machines may not be able to obtain the same results as the conventional method. Therefore, we report a new method of the phase-field model whose final simulation state must be the global minimum free energy.
To perform phase-field simulation by the Ising machine, we transform the phase-field model into a QUBO model. We set a two-dimensional calculation area and divide it into grids. We assign four binary variables to any i-th grid point S i : and we express the order parameter of the i-th grid point as 0≤k<4 a 4i+k . Since a 4i+k ∈ {0, 1} , the value range of S i is from 0 to 4 in this case. If we assign more binary variables to each grid, the resolution of the result is expected to increase, and as a result, the total number of required binary variables also increases.
First, we conduct numerical simulation cases with and without Ohta-Kawasaki energy. We set the material properties in the analysis as follows: α I = 20.0, α A = 5.0, α OK = 0.0/0.4 and f = 0.5 . Figure 3 shows the numerical analysis results of S i . Clearly, the phases of polymer A and polymer B, shown in white and black, respectively, www.nature.com/scientificreports/ are separated into two large phases without Ohta-Kawasaki energy, as shown in Fig. 3a, while the composition of polymers A and B with Ohta-Kawasaki energy equilibrates with a finer pattern, as shown in Fig. 3b. Figure 4 shows the emulation process of the annealing optimization in other cases with and without Ohta-Kawasaki energy. Note that this figure shows the results obtained when the Ising machine is interrupted before the process has finished and does not directly correspond to the physical time in the conventional method. We set α I = 20.0, α A = 5.0, α OK = 0.0/0.5 and f = 0.3. The simulation result without Ohta-Kawasaki energy (Fig. 4a(iv)) shows a larger pattern than the one with Ohta-Kawasaki energy (Fig. 4b(iv)). With the proposed Ising models, the simulation takes only 1 s.
To show the exact execution times of the two methods, we studied how the energy decreases with time in the conventional method and our new method in one simulation as an example. In Fig. 5a, we show the time evolution of the free energy of one randomly selected simulation using the conventional method. This figure shows that the simulation using the conventional method converges in about 4000 s (64 × 64 cells). Correspondingly, , the observed microstructure corresponds to that in the phase diagram, i.e., a lamellar structure. www.nature.com/scientificreports/ we performed some simulations using our new method with different timeout times. In Fig. 5b-d, we show the minimum energy that was found in less than a specified time. This figure shows that the simulations using Ising machines reach the minimum energy in about 0.8 s (32 × 32 cells), 1.5 s (40 × 40 cells), and 2.5 s (48 × 48 cells), which is very fast compared with the conventional method. Finally, we carry out inclusively large-scale parametric analyses. Figure 6a shows the numerical analysis results in the case of fraction rate f = 0.5 , and Fig. 6b shows the case of f = 0.3.  www.nature.com/scientificreports/ When you look at the phase diagram in Fig. 1, the case of f = 0.5 shows only lamellar patterns, indicating that the parametric study using the defined Hamiltonian outputs physically correct results. In the case of f = 0.3 , hexagonal or cubic patterns should be observed from the phase diagram of Fig. 1, and it seems that the dotted pattern or the striped patterns may change depending on the viewing direction of these patterns because of the carrying out of two-dimensional calculations. Similar to the conventional phase-field analysis results, the characteristic length of the pattern seems to be shortened; thus, it makes sense that the pattern becomes finer as the long-distance interaction becomes stronger with increasing Ohta-Kawasaki energy.
When the coefficient of the gradient energy term becomes large, the parameter corresponding to the diffusion coefficient increases. This parameter change does not appear in the phase diagram in Fig. 1. The phase diagram should generally be created with the same material at a constant environmental temperature. However, the above outcome makes sense because the larger the gradient energy becomes, the larger the structure tends to be (the farther the material is distributed) in Fig. 6a,b. Additionally, in the region where the interaction energy is small, the pattern tends to be disordered, which is consistent with the phase diagram.
As shown in Fig. 6, the results are quite similar to those obtained by the conventional method, which means that the proposed modeling based on the Ising model is quite useful.
Note that there is no strict symmetry in the Ohta-Kawasaki term, but essentially, in the case of f = 0.7 and f = 0.3 , the colors are simply swapped. If α OK equals zero, there is strict symmetry; thus, the same pattern with completely inverted colors will appear. Additionally, since the Hamiltonian has symmetry when rotated 90 degrees about the xy axes, the direction of the lamellar structure appears randomly.
In terms of accuracy evaluation, the new QUBO model is not exactly the same as the conventional one and, thus, does not shared exactly the same parameters as the conventional method; making a direct comparison with the same model parameters is impossible. Instead of comparing the accuracy between the two methods, we evaluate errors due to the discretization of continuous order parameters in a simulation. In our QUBO model, the order parameters S i was discretized to (0, 1, 2, 3, 4). We now introduce a continuous version of our QUBO model, called continuous-QUBO, that is, the term k a 4i+k in Eqs. (11)-(14) is replaced with continuous variables x i ∈ [0, 4] . Then, the Hamiltonian of continuous-QUBO is minimized and the optimal order parameter x * i is obtained using the gradient descent method. The error is calculated between the result obtained with the Ising machine, S * i , and the result of continuous-QUBO, x * i , which is considered the true value. This enables us to calculate the accuracy of the same model parameters. Figure 7 shows (a) one of the results obtained with the Ising machine, S * i , (b) the corresponding optimal order parameter of continuous-QUBO, x * i , and (c) the difference between the two plots. In this result, the mean error per cell � S * i − x * i � i is about 0.018, which is sufficiently smaller than 1.0, which is the discretization size of the order parameter S * i . In this study, we proposed a novel method to solve the phase-field model by the objective function for exploring the phase-separated structure of diblock polymers by an Ising machine, with an eye toward the use of quantum annealer. We confirmed the tendency from the obtained results and the performance by comparing them with the phase-field simulation. The approach involves the use of a global optimization metaheuristic algorithm, called simulated annealing, to directly minimize the Helmholtz free energy instead of minimizing it analytically and then solving the resulting nonlinear, partial differential equation, i.e., the governing equation of the phase-field model. As a result, we obtained the following: 1. The results obtained by Ising machine showed the same tendency as the results obtained from the conventional phase-field analysis. 2. The obtained results were consistent with the phase diagram. 3. The analysis time was shortened to the extent that high-speed comprehensive analysis becomes possible.
The proposed modeling based on the Ising model has been shown to be very useful. In the future, we will address the following remaining issues for the practical use of this method. First, we will investigate the usefulness of the intermediate results obtained with an Ising machine, which may correspond to the non-equilibrium a b c www.nature.com/scientificreports/ time history of the distribution obtained by the conventional method. Second, as the number of bits that Ising machines can handle increases, it will be possible to apply our method to the 3D domain. Third, as noted above, a direct comparison with the same model parameters is not yet possible. We will clarify the relationship between the parameters of our model and these of the conventional model by introducing correction terms using machine learning methods. This scheme on how to discretize the continuum variable and formulate the QUBO will help in solving the other problems related to continuum mechanics.

Methods
Theory of the phase-field model. The phase-field model for spinodal decomposition is introduced. First, the Cahn-Hilliard equation 31 can be written as follows: where c is the concentration of a phase, which is the order parameter in the phase-field model, F is the free energy functional and M is the mobility. Here, F is formulated with consideration of the long-distance interaction as follows: where F grad is the gradient energy, F chem is the chemical potential and F long is the long-distance interaction. Here, µ is defined as µ ≡ −κ∇ 2 c + χ RTlnc , and M is considered constant. For the implementation of Eq. (8) in a program code, the finite difference method is generally employed. In addition, , and δ is the Green function, which is calculated through a Fourier transformation.
Calculation by an Ising machine. We conducted the calculations using Fixstars Amplify Annealing Engine (Amplify AE) 32 with timeout of 1 s. Amplify AE is GPU-based Ising machine that can handle 100,000 bit-class problems.
When solving a problem by using Ising machines, we must design QUBO models of the binary variables. The QUBO model is formulated as Hamiltonian. We formulate the Hamiltonian. We construct the whole Hamiltonian H as a linear combination of the summation preservation term H sum , the interaction (internal) energy term H int , the gradient (adjacent) energy term H adj , and the long-distance energy term (the so-called Ohta-Kawasaki energy) H long : where The term H sum constrains the total order parameter of space. Equation (11) uses the square constraint because a strict preservation term cannot be used in QUBO models. The constant α F controls the strictness of preservation, and the constant f denotes the target total ratio of the order parameter. The term H int makes the order parameter of each point 0 or 4 to avoid intermediate states as much as possible, with strength α I . The term H adj controls the strength of adjacent interactions by bringing closer the order parameter of adjacent pairs of the grid points together. P adj denotes the set of all adjacent pairs. The coefficient α A is the strength of adjacent interactions. The term H long expresses the Ohta-Kawasaki energy, and g ij is defined as g ij ≡ 1/r ij , with the distance between grids i and j denoted by r ij . P all denotes the set of all pairs of all grid points.
Multiplying the coefficient of each term by a constant does not influence the simulation results; thus, we set α F to 1. In this setting, the whole Hamiltonian contains at most a quadratic term of the variables {a i } ; therefore, rearranging the whole Hamiltonian immediately yields the coefficient of QUBO models Q ij .

Data availability
The datasets used and analysed during the current study available from the corresponding author on reasonable request.